数学模型

Wang Haihua

🍈 🍉🍊 🍋 🍌

Logistic回归模型

模型的构建

Logistic 回归, 与普通回归任务不同, 分类的结果输出是有限的离散值。 以二分类为例, 结果输出要么为 0 , 要么为 1 , 即 $y \in\{0,1\}$ 。其基本思想是 在空间中构造一个合理的超平面, 把空间区域划分为两个子空间, 每一种类 别数据都在平面的某一侧。不能按照多元线性回驲式 $$ f(x)=\beta_{0}+\beta_{1} x_{1}+\cdots+\beta_{m} x_{m} \text { (12.18) } $$ 来分类。

因为直接使用式 (12.18) 得到的是实数值, 需要将实数值规约为 0 或 1, 较理想的是阶跃函数 $$ f(x)= \begin{cases}0, & \beta_{0}+\beta_{1} x_{1}+\cdots+\beta_{m} x_{m} \leq 0, \\ 1, & \beta_{0}+\beta_{1} x_{1}+\cdots+\beta_{m} x_{m}>0 .\end{cases} $$ 但阶跃函数在 0 点不连续、不可导。为此, 可以通过阶跃函数的平滑版本, 即 sigmoid 函数来为我们实现: $$ f_{\beta}(x)=\frac{1}{1+e^{-\left(\beta_{0}+\sum_{j=1}^{m} \beta_{j} x_{j}\right)}}, $$ 其中 $\beta=\left[\beta_{0}, \beta_{1}, \cdots, \beta_{m}\right]^{T}$ 。

sigmoid 函数具有很多优秀的性质：首先, 它将输入数据压缩为 0 到 1 的范围内, 得到的结果不是二值输出, 而是一个概率值, 通过这个数值, 可 以确定输入数据分别属于 0 类或属于 1 类的概率: $$ \begin{gathered} P\{y=1 \mid x\}=f_{\beta}(x)=\frac{e^{\beta_{0}+\sum_{j=1}^{m} \beta_{j} x_{j}}}{1+e^{\beta_{0}+\sum_{j=1}^{m} \beta_{j} x_{j}}}=p \\ P\{y=0 \mid x\}=1-f_{\beta}(x)=\frac{1}{1+e^{\beta_{0}+\sum_{j=1}^{m} \beta_{j} x_{j}}}=1-p . \end{gathered} $$

由式 (12.21) 和 (12.22) 可以得到 $$ \ln \frac{p}{1-p}=\ln \frac{P\{y=1 \mid x\}}{P\{y=0 \mid x\}}=\beta_{0}+\sum_{j=1}^{m} \beta_{j} x_{j} . $$ 式 (12.23) 称为 Logistic 回归模型, 式 (12.23) 的左边是与概率相关的对 数值。所以, 无法使用通常的最小二乘法拟合末知参数向量 $\boldsymbol{\beta}$, 而是采用极 大似然估计法。 式（12.23）也可以改写为 $$ p=\frac{1}{1+e^{-\left(\beta_{0}+\sum_{j=1}^{m} \beta_{j} x_{j}\right)}}, $$ 其中, $p$ 为事件发生的概率。

Logistic 模型的参数估计

为了构造似然函数, 我们把式 (12.21) 和 (12.22) 统一改写为 $$ P\{y \mid x ; \beta\}=f_{\beta}(x)^{y}\left(1-f_{\beta}(x)\right)^{1-y}, y=0 \text { 或 } 1 . $$ 为了拟合模型中的末知参数向量 $\boldsymbol{\beta}$, 使用极大似然估计法, 需要构造似 然函数。似然函数的统计背景是, 如果数据集 $\left\{\left(x_{i}, y_{i}\right), i=1,2, \cdots, n\right\}$ 中每个样 本点都是相互独立的, 则 $n$ 个样本点发生的联合概率就是各样本点事件发生 的概率乘积, 故似然函数可以表示为 $$ L(\beta)=\prod_{i=1}^{n} P\left\{y_{i} \mid x_{i} ; \beta\right\}=\prod_{i=1}^{n}\left[f_{\beta}\left(x_{i}\right)^{y_{i}}\left(1-f_{\beta}\left(x_{i}\right)\right)^{1-y_{i}}\right] . $$

为了求解方便, 将似然函数做对数处理, 得到 $$ \tilde{L}(\beta)=\sum_{i=1}^{n}\left[y_{i} \ln f_{\beta}\left(x_{i}\right)+\left(1-y_{i}\right) \ln \left(1-f_{\beta}\left(x_{i}\right)\right)\right] . $$ 拟合参数向量 $\boldsymbol{\beta}$, 转化求函数 $\tilde{L}(\boldsymbol{\beta})$ 的最大值, 无法求得其解析解, 只能求数 值解, 我们可以使用经典的梯度下降算法求解参数向量 $\boldsymbol{\beta}$, 具体算法我们这 里就不给出了。下面给出在分组数据情形下参数向量 $\boldsymbol{\beta}$ 的最小二乘估计。

分组数据情形下参数的最小二乘估计

在对因变量进行的 $n$ 次观测 $y_{i}, i=1,2, \cdots, n$ 中, 如果在相同的 $x_{(i)}=\left[x_{i 1}, x_{i 2}, \cdots, x_{i m}\right]^{T}$ 处进行了多次重复观测, 这种结构的数据称为分组数 据, 分组个数记为 $c$; 则在式 (12.23) 中可用样本比例对概率 $p$ 进行估计, 对应 $x_{(i)}$ 的概率估计值记作 $\hat{p}_{i}$, 并记 $$ y_{i}^{*}=\ln \left(\frac{\hat{p}_{i}}{1-\hat{p}_{i}}\right), i=1,2, \cdots, c, $$ 则得 $$ y_{i}^{*}=\beta_{0}+\sum_{j=1}^{p} \beta_{j} x_{i j}, i=1,2, \cdots, c . $$

由线性回归模型的知识知, 参数 $\beta=\left[\beta_{0}, \beta_{1}, \cdots, \beta_{m}\right]^{T}$ 的最小二乘估计为 $$ \hat{\boldsymbol{\beta}}=\left(X^{T} X\right)^{-1} X^{T} y^{*}, $$ 其中 $$ \boldsymbol{y}^{*}=\left[\begin{array}{c} \boldsymbol{y}_{1}^{*} \\ \boldsymbol{y}_{2}^{*} \\ \vdots \\ y_{c}^{*} \end{array}\right], X=\left[\begin{array}{cccc} 1 & x_{11} & \cdots & x_{1 m} \\ 1 & x_{21} & \cdots & x_{2 m} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{c 1} & \cdots & x_{c m} \end{array}\right] . $$ 下面用一个例子来说明分组数据 Logistic 回归模型的参数估计。

案例

问题

在一次住房展销会上,与房地产商签订初步购房意向书的共 有 $n=325$ 名顾客，在随后的 3 个月的时间内，只有一部分顾客确实购买了 房屋。购买了房屋的顾客记为 1 , 没有购买房屋的顾客记为 0。以顾客的家 庭年收入为自变量 $x$, 家庭年收入按照高低不同分成了 9 组, 数据列在表 $12.4$ 中。表 $12.4$ 还列出了在每个不同的家庭年收入组中签订意向书的人数 $n_{i}$ 和相应的实际购房人数 $m_{i}$ 。 房地产商希望能建立签订意向的顾客最终真正 买房的概率与家庭年收入间的关系式，以便能分析家庭年收入的不同对最 终购买住房的影响。

\begin{aligned} &\text { 表 } 12.4 \text { 签订购房意向和最终买房的客户数据 }\\ &\begin{array}{cccccc} \hline & \text { 家庭年收入号 } & \text { 签订意向书 } & \text { 实际购房 } & \text { 实际购房比例 } & \text { 逻辑变换 } \\ & \text { (万元） } & \text { 人数 } n_{i} & \text { 人数 } m_{i} & \hat{p}_{i}=\frac{m_{i}}{n_{i}} & y_{i}^{*}=\ln \left(\frac{\hat{p}_{i}}{1-\hat{p}_{i}}\right) \\ \hline 1 & 1.5 & 25 & 8 & 0.32 & -0.7538 \\ 2 & 2.5 & 32 & 13 & 0.4063 & -0.3795 \\ 3 & 3.5 & 58 & 26 & 0.4483 & -0.2076 \\ 4 & 4.5 & 52 & 22 & 0.4231 & -0.3102 \\ 5 & 5.5 & 43 & 20 & 0.4651 & -0.1398 \\ 6 & 6.5 & 39 & 22 & 0.5641 & 0.2578 \\ 7 & 7.5 & 28 & 16 & 0.5714 & 0.2877 \\ 8 & 8.5 & 21 & 12 & 0.5714 & 0.2877 \\ 9 & 9.5 & 15 & 10 & 0.6667 & 0.6931 \\ \hline \end{array} \end{aligned}

解 显然, 这里的因变量是 0-1型的贝努利随机变量, 因此可通过 Logistic 回归来建立签订意向的顾客最终真正买房的概率与家庭年收入之 间的关系。由于表 $12.4$ 中，对应同一个家庭年收入组有多个重复观测值, 因此可用样本比例来估计第 $i$ 个家庭年收入组中客户最终购买住房的概率 $p_{i}$ 其估计值记为 $\hat{p}_{i}$ 。然后, 对 $\hat{p}_{i}$ 进行逻辑变换。 $\hat{p}_{i}$ 的值及其经逻辑变换后的值 $y_{i}^{*}$ 都列在表 $12.4$ 中。

本例中, 自变量个数 $m=1$, 分组数 $c=9$, 由 (12.18) 式计算可得 $\beta_{0}, \beta_{1}$ 的最小二乘估计分别为 $$ \hat{\beta}_{0}=-0.8863, \hat{\beta}_{1}=0.1558, $$ 相应的线性回归方程为 $$ \hat{y}^{*}=-0.8863+0.1558 x, $$ 决定系数 $R^{2}=0.924, F$ 统计量 $=85.42$, 显著性检验 $p$ 值 $\approx 0$, 线性回归方程 高度显著。最终所得的 Logistic 回归方程为 $$ \hat{p}=\frac{1}{1+e^{0.8863-0.1558 x}} . $$

由 (12.29) 式可知, $x$ 越大, 即家庭年收入越高, $\hat{p}$ 就越大, 即签订意 向后真正买房的概率就越大。对于一个家庭年收入为 9 万元的客户, 将 $x=x_{0}=9$ 代入回归方程 (12.29) 中, 即可得其签订意向后真正买房的概率 $$ \hat{\boldsymbol{p}}_{0}=\frac{1}{1+e^{0.8863-0.1558 x_{0}}}=\mathbf{0 . 6 2 6 2} \text {. } $$ 这也可以说, 约有 $62.62 \%$ 的家庭年收入为 9 万元的客户, 其签订意向后会 真正买房。 把表 $12.4$ 中 $x, n_{i}, m_{i}$ 三列数据保存到文本文件 Pdata12_7.txt 中。

代码

import numpy as np
import statsmodels.api as sm
a=np.loadtxt("data/house.txt")   #加载表中x,ni,mi的9行3列数据
x=a[:,0]; pi=a[:,2]/a[:,1]
X=sm.add_constant(x); yi=np.log(pi/(1-pi))
md=sm.OLS(yi,X).fit()  #构建并拟合模型
print(md.summary())  #输出模型的所有结果
b=md.params  #提出所有的回归系数

p0=1/(1+np.exp(-np.dot(b,[1,9])))
print("所求概率p0=%.4f"%p0)
np.savetxt("data/house2", b)  #把回归系数保存到文本文件